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This  report  presents  a  sunmary  of  the  aathenatleal  equations  contained 
in  the  R-SIMS  missile  simulation  program.  This  program^  vrltten  for  the 
IBM  T090  digital  c(miiputer>  is  designed  to  he  a  flexible,  hl|^  speed  computation 
tool,  useful  for  guidance  system  analysis,  trajectory  design,  hardvare 
evaluation,  or  post-flight  enalysls  of  any  current  or  proposed  missile  system. 

The  simulation  is  generally  complete  and  in  most  instances  accuracy  vlU  he 
limited  only  hy  the  accxiracy  of  the  input  data. 

Currently,  vork  is  proceeding  vhlch  vill  add  the  effects  of  missile 
rotational  inertia  to  the  simulation. 

The  mathematical  model  'ised  is  fundamentally  simple.  The  forces  acting 
on  an  assumed  "point  mass"  are  computed  and  divided  hy  the  mass  of  the  missile 
to  obtain  accelerations.  These  accelerations  are  added  vectorially  to  that  of 
gravity  and  numerically  integrated  to  obtain  a  velocity  increment  over  a  short 
period  of  time.  The  velocity  is  numerically  integrated  to  obtain  a  position 
increment  for  the  same  period  of  time.  As  this  process  is  repeated,  the  normeJ. 
output  of  the  simulation  is  the  tine  history  of  a  trajectory. 

Di:icor.tinultie8  in  the  trajectory,  such  os  stagings,  Jettisons,  etc.,  are 
functions  of  input  data  and  are  treated  logically  by  the  program. 

In  addition  to  the  quantities  necessary  to  the  simulation,  other  descriptive 
quantities  are  computed  and  printed. 

For  design  and  post-flight  analysis  problems,  automatic  iterations  (over 
traject’  rle*:  jr  part.  1  trajectories)  may  be  made  by  the  program.  The  iteration 
techniques  aliev  any  nunber  cf  input  parameters  to  be  found,  with  reference  to 
values  for  any  number  of  fundamental  or  descriptive  quantities.  Maximization  or 
minimization  may  also  be  accomplished,  with  or  without  the  presence  of  other 
constraints. 

Although  numeric ed  integration  is  at  the  heart  of  the  program,  no  statement 
of  the  Integration  techniques  used  is  here  made.  The  Interested  reader  is  referred 
to  any  standard  numerical  analysis  text  for  a  discussion  of  Runge-Kutta  and 
Adams -Moulton  Integr^ition  techniques. 
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Derivatives  of  Equations  of  Motion 

All  of  the  basic  confutations  are  carried  out  In  a  Cartesian  Inertial  system 
having  its  origin  at  the  center  of  the  reference  ellipsoid.  The  x  and  y  axis  pass 
through  the  equator  and  the  z  axis  is  along  the  polar  axis.  (Fig.  l) 

Msslle  attitude  determines  the  directions  of  the  force  vectors.  Mutually 
perpendicular  unit  vectors  are  defined  such  that  £  points  along  the  missile  roll 
axis,  7^  lies  in  the  pitch  plane,  and  ^  lies  in  the  yaw  plane.  (Fig.  l) 

Thrust  and  Moss  Flow 
sure  given  by 

where 

F^,  may  result  from  a  table  interrogation  or  from  the  Atlas  Influence 
Coefficient  Engine  Model. 

Velocity  relative  to  air  mass  is 


where  wind  velocity  components  are 


Derivatives  of  Equations  of  Motion 


and 


-V 

w 


K 


cosA 

z 


(north  cojnponent  of  vind) 


W  *  -V  K  sinA  (east  component  of  wind) 

0  V  p  2* 

Atmospheric  quantities  p  (alugs/ft'^),  C  (ft/sec),  T  (deg  Rankine)  and  P  (Ihs/in  ) 
are  determined  hy  approximation  to  the  ARDC  model  of  1959  ®nd  are  avallahle  to  2.3  million 
feet  of  altitude. 

Radius  at  sea  level 


where 


. 


.  2 
K 


2e  - 
(1  -  e) 


ellipticity  (flattening) 


AM  li'idr 


-  ^SL 
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Derivatlvea  of  Equations  of  Motion 


Thi-^ist  and  Mass  Flov  (Cont'd) 


Mtich  niuater 


Dynamic  pressure 


Aerodynamic  Drag 


=  S  «  s 

where  is  obtained  from  a  table,  normally  as  a  function  of 
Center  of  Gravity  Offset  Force 

(See  Fig.  3)  is  obtained  directly  from  a  table  or  i®  obtained  from  a 

table  and 

^CG  "  ^  ■  ^CG 

where  is  nose  to  center  of  gravity  distance  along  ^  axis ■  The  CG  distances 

in  the  direction  and  ^  direction  (^^q)  are  tabular. 

Center  of  gravity  offset  forces  due  to  drag  in  7J  and  ^  directions  ere 


^'£0  'I 


(-7-) 
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Derivatives  of  Sq:uatlons  of  Motion 


Center  of  Gravity  Offset  Force  (Cont’d) 


Center  of  gravity  offset  forces  due  to  thrust  nlsaUgODent  In  7|  and  ^  directions 


P  >  "  A  ^"7 — 

^  P  h  F  ^CG 


Aerodynamic  Normal  Force 


Pitch  and  yaw  components  of  angle  of  attack  (a),  angle  "between  £  and  aie 


determined  from 


tan  a  =  - 
P 


r}'\  7? .  V 

-=--3r  =  ' 

£•  \  \\\ 


(the  approximate  form  is  useful  for  mnan  values  of  o) 


tan  a  =  - 

y 


~  *  2 
a  =  tan  a 
p 


+  tan 


If  the  approximation  is  used,  o:  and  or  are  assumed  equal,  respectively,  to  tan  ct  and 

P  y  p 


t  or.  a  , 


7]  (yav  axis) 


FIGURE  2 
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Derivatives  of  Iquatloos  of  Notion 


Aerodynamic  Nonnal  Force  (Cont'd) 

Aerodynamic  normal  force  coefficients  in  the  pitch  and  yav  planes  (C„  i  C„  ) 

may  he  pxresented  as  tabular  functions  of  a  .  a  and  N  .  When  C„  and 

y  m  N 

C  are  presented  ® 


S  ■  S  “p 

p  Po 


y 


or,  a  tahxilar  normal  force  coefficient  may  he  entered  directly  or  as  a  function  of 
total  angle  of  atteuik: 

\  “ 

tan  a 

S  “  tan 
P 

tan  a 
S  "  tan 

y 

Either  £  are  obtained  directly  from  a  table  or  £  are  obtained 


from  a  table  and 


^  --  L  -  S 

'  cp  ,  -^cp  I 

ar:'  ^he  nose  to  center  of  pressure  distances  in  the  pitch  and  yav  planes 


^  a:-;l  ■  (Pip^. 
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Dexivatlvet  of  Equations  of  Notion 

Gravity  (Coot'd) 


and  J,  E  and  P  are  constants,  and  optionally  need  not  l>e  recomputed  for 
aoall  changes  in  z  position. 


Acceleration 


in  the  case  of  spherical  body  drag.  In  this  case  ablation  effects  may  be 

{ 

Included i 

"-"I'm 
®  ■*!  *8 

and  Rg  are  tabular  mass  and  cross-sectional  area  ratio  factors. 


vhere  R^ 
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DerlTBtlTes  of  SguctlaoB  of  liotlaD 
Integral*  for  Velocity,  Position  and  Mmb 

'v  =  J  V  dt 

't  ^  y  V  dt 

M  =  J  M  dt 

Turning  Rate  Options 

Several  computed  turning  rate  options  are  avallaible.  In  each  option  the  rate 
computed  (os  )  is  limited  to  +  2°  per  integration  cycle. 

Zero  lift  flight  (gravity  tvm) 

Initially  orient  V  to  obtain  a  *  0 


V  =  V 

a  a 


l! 


then 


Q 


+  V  +  V 
a  V 


end  comrute  pitch  and  yaw  rates  to  med-ntain  a  constant  attitude  relationship 

_Jb 

between  £  and 


OJ 


Pg 


yg 


ri-v. 

_ I  ~ 
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Do^v»tlves  of  BquAticxQB  of  Motion 


(Cottt»d) 


irtiero 


(x  +  y  ) 


.)  *  <V>‘ 


VW2  W  ^ 

- - ^T72  +  2  Q  V  +  (Q  X 

2  2//2  e  vy  '•’‘e'' 

(jc  +  y  ) 


VWl  w 


(x"  ^  y^) 


TI^  -  V„  *  (Q^)2  y 


rrry 


(Wy)  X  z 


1"?!  ‘  ( 


2  2. 
X  +  r  ji 


(H^)  y  1 

iTTl  i? 

I  r[  (x  +  y  ) 
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Derlvfttivee  ot  Icufttii 


Turning  Rate  Options  (Cont'd) 

Roll  to  Qedntaln  pitch  plane  through  center  of  Cartesian  Inertial  ajreteni 


U) 


rg 


■TT 


U) 


note  that  is  indetenninant  during  vertical  rise. 
Instantaneous  turns  about  pitchy  yaw  and/or  roll  axes  (Fig.  2). 

A  =  cosO  -  ^  sin© 
b  o  p  'o  p 

7?  =  £  sin©^  +7?  cos© 

(  o  p  '  o  p 


^  =  £  cos©  +  £  sin© 
o  o  y  o  y 


t  ^  -  C  sin©  +  C  cos© 
3  o  y  o  y 

?“  £ 

V  n 


7?=  7?^  cos©^  .  ^ 


sin© 
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Derivatives  of  Equations  of  Motion 


Tumlnp  Rate  Options  {Cv..it'd) 

Obtain  a  specified  epsilon  local  (Cj.  )(see  belov)  in  a  specified  time  (At) 
from  the  current  ) 


where 


€j^  is  the  euigle  between 


£  and  r. 


(Fig.  6) 


Obtain  €,  imnedlately  from  €_ 
*13  c 

I-  ?  -  n 

JO  o 

n-  X  »in4' .  ccf 


Make  €_  =  eaigle  between  r  and  V  (Fig.  5) 
Xi  X 


■  % 

“pg  “  6t  ^  *  'l 


where  Ot  is  the  computing  Incresaent  and 


|v|  SU,  gj 

r 
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vbere 
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Auxillftry  CoBiputfttloiJO 

Descriptive  Quantities 

These  e<).vistloQ8  specify  some  descriptive  guantltles  used  for  special  puiposes 
aiid/or  print  cmt. 

Angle  hetween  V  and  'r  (Fig.  4) 

& 


p  B  COS 


-1  \ 


It! 


0  *  p  2  i8o 


r- 1  -  ^ 

Angle  hetveen  ^  end  V  (Fig.  3) 


Pj  «  cos 


-1  r  ♦  r 
-* 

r  r 


0  i  2  IBO 


ri'l-pi 


Angle  hetveen  "r  and  £  (Fig.  6) 


=  cos 


0  2  €j^2  180 


r  •  ^ 

.Angle  hetveen  ^  and  launch  vertical.  (i  )  (Fig.  7) 


C  B  COS 


0  2  6  2  100 


Azimuth  from  north  of  S  projected  on  a  plane  normal  to  r 


^body 


tan 


c  p 
X  -  y 


r— 

^  p - T 

y  S  -  z  $ 

-  X 

z  S .  -  X  S„ 

^  z  y 

X  z 

Auxiliary  ConnputationB 
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blve  Quantities  (Cont'd) 

Nonoallzed  angular  noisentuoi 
RVSB  -  |t|  |v|  Bin  Pj 

Ratio  squared  of  velocity  to  circular  orbital  velocity 


To  tal  energy 


Vis  viva  ei»ergy 
WEN  •  2  KNGY 
Sezised  velocity 

Vg  =■  /  (V  -  G)  dt 
Aerodynamic  heating 

Geodetic  latitude  (e^)  and  longitude  (J^)(Fig.  8) 


X 

lATrrUDE  i  IDNGITUDE 
FIGURE  8 
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Auxiliary  Confutations 

I 

Descriptive  Quantities  (Cont'd)  ;  Dei 


o  s  tan 
c 


-1 


L (x^  +  y2) 


Circular  range  (r)  tetween  two  positions  on  the  spheroid^  where  (9^^, 
(geocentric  latitude,  longitude)  defines  the  first  position, (0^2 ,  defines  the 
second,  and  r^  is  the  radius  of  the  spheroid,  is  given  by 


If 


=  ®^®cl  ‘=°®®cl  ‘'°®®c2  -  ^1^ 

2:  .99998,  then 


where 


If 


R  = 


.99998,  then 


2x  -  y,  if  90°  IS 

^  ,lfO°=E 


180° 

90° 


Htsfi 


I  tota 


where  p  is  the  clockwise  difference  of  the  azimuth  frciB(0  ,,  ^.)  to  (O  «,  and 

CX  X  CC  b 


«  reference  azimuth. 
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Auxiliary  CoiqputatloQs 


Mse  Dlstaacea  (Coot'd) 

If 

“l,I  “  -^9998,  than 

*  [<®CI  -  *  <■*!  ■  ^1?“ 

|®CI  -*Cl|“*|^I  "*  =  ^° 

If 

•*  *99996,  then 

T  *  cob'^JS-  _:  0  —  ^  S IfiO^ 
L,I  ■  *L,I"  ^  L,I 


\h  *  ®^®CL  ®^®CM 
If 

®L,M  “  -99998,  then 
\m  =  [^®CM  -  *CL^^ 

were 

|®CM  *  ®Cl  «4  I  ^  [  ««  =  180° 

If 

®L,M  -99996,  then 


■  COB 


2180° 
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Auxiliary  Conputa-tlons 


Miss  Distances  (Cont’d) 


M.  =  m; 

c  o  c 


«t  '  ^  "d'’]  ’ 

vhere  M',  M'  are  central  angles.  Dote  that  la  poaltive  to  tbe  left  vhen 

c  d 

looking  dovn  range. 

Orbital  and  Impact  Approximations 

For  on  elliptical  orbit  or  an  impact  trajectory, 


^  m  •  ivr  T 

c  cosE  =  * — '  p  “ 


1  .  e  cosE 


e  sinE  = 


r  •  V 
A  •  G,, 


(e  cosE)^  +  (e  sinE) 


Perigee  =  A(1  -  ^)"^q 


Apogee  --  A(1  + 


£  Is  the  ellipse 


idol  eccentricity  and  E  is  the  eccentric  anomaly  of  the 


.  ..  g. ■  i i', fj  a.  1  c  posi.*  ion. 


AxixlXlary  Conmutatlax^ 


Orbital  «nd  Impact  Approxlawtlcaae  (Cont’d) 


Period  ■ 


2,  A  3/^ 


Inclination  ■  cos"^  (slxiA^  coaO^j^) 


Velocity  at  apogee 


1^1  |v|  (1  -  cos^pj) 


1/2 


A(1  +  c) 


For  an  isqpact  trajectory  Rgg^  2:  A(l  -  e) 


c  coal^  •  1 


aei 


r  2  2I 

c  ain^j  -  -  j^€  -  (€  cosl^)  J 


1/2 


(c  sln^)(c  cosE)  -  (c  coB^j,)(c  sin^ 

Bin(]^  -  K)  -  - i- - g- - T - 


_  (a  coa]^)(fi  cosi)  +  («  8i]3]l^)(c  slnf) 
C08(^j  -  *)  •  - i— - -g — : - i-T - — — 


f  *  tan 


-1 


8ln(l^  -  I) 

(1^  -i) 


COS 


xS  wiO  €CC£ii‘bx^j.C  iujOvniU^  XQ^ttCC^  pOXXl*C* 


Page  27 


Auxiliary  Ccuqputatloas 


Lons  (Cont'd) 


Ti»e  of  flight  to  impact  and  time  of  iu^pact  are 


k 


S  +  c  sinE  -  €  2 


T  =  t  +  T 


Inertial  Cartesian  con^nents  of  inspect  point  are 


*aEi  ■  *•  *8 


^  *  >8 


'SK  ■  *  "8 


vhere 


r  3/2  r  ,  ~  1 

I  -  B)  ♦  *  sin*  -  * 

co8(£,  -  E)-e  cosE 

g  =  - - - 

1  .  e  cosE 


Azimuth  of  impact  Is 


-  ^tan”^ 


8in(l80  -  A,  ) 


-  iCO 
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Auxllisory  Goarpu't-s'tioiis 

Orbital  and  lapoct  Approximations  (Cont’d) 
where 

cos©  6ln(0  -  ^.) 

sin  (180  -  A.  )  =  - Si - ^ - S- 

elnY 


ein  e  .  cobH^  -  sin  0 

cos  (i8o  -  A  )  - - nr"— ^ 

cos  alnT 


If 


«  slnd^j.  8ln0^,j  +  cose^,^  cosO^^  cob(;^j^  -  0^) 


=  .99998,  then'l'  .  [(8„  -  i  (|<j  -  oos\^_ 


where 


|«ci  -®cl!  |»‘i  -’'lI  »«=18o' 


If 


-c  .99998,,  ^ 

Circular  range  to  impact  is 
2k  “ 

CRIP  = 


180 


j  if  90^  2 

jp|  =  180' 

;  if  0°  2 

|p|  -  90' 

1/2 


where  p  is  the  clockwise  difference  of  ii^pact  azljoutb  and  launch  azimuth. 
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Auxiliary  Ccsput^ti 

Orbital  and  Intpact  /feproxlroattona  (Cont*d) 

If 


I”?  ^  or  if  the  Tile'',i\e  is  in  a  hyperbolic  oxhit  {i.e. 


MJlli.j 


'M 


2),  all  orhital  od  impact  quantitiaa  are  aet  to  »ero. 


Rggj  (radius  of  spherical  earth  impact)  may  be  input  directly  or  as  a 
hel^t  relative  to  the  radius  of  the  target  ^Rggj  ■  or  as  a 

hel^t  rel  ative  to  the  everage  radius  of  the  geoid  (r^)» 

Velocity  Losses  and  Impulse  Calcul-ationg 

Gravity  loss 


f 

£ 


r 

=■  I  g  cos^j  dt 


vhere 


=  time  at  beginning  of  stage 
tg  *  time  at  end  of  stage 
g  =  local  gravity 

-A  ^ 

Pj.  =  angle  between  r  and  geocentric  local  vertical  (1) 


Drag  loss 


AuxtUtey  Cc!?S!*it-ti  AVW* 

Velocity  toaecB  and  InsnOiie  Calculatlona  (Coot’d) 
Arerege  ■peciflc  tiqpulfle  j)er  stage 


Total  iayulae 


Ideal  velocity  per  stage 
•  ®o  ^ap 

vbere  ^  1.  the  nominal  vnlun  of  grmvlty,  i,  tb«  wtgfct  of  mloslle 

at  and  vtg  la  the  weight  of  mlaalle  at  t  . 

6 
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t 


i 


F«^  31 


Distance  from  radar  site  to  polar  axis  la 


■^SI^  ’^at  *  <“1!  *  «aR>  “•^.J 


vhere 


X  "  1 - 1 — r 

^  1  +  (l-e)'  tan 


a  altitude  of  radar  site  aibove  sea  level 


at  geoidal  8epar«itioa  at  radar  site 

Tor  *  level  radius  at  radar  site 

00-  ■  longitude  of  radar^  positive  east 
li 

Distance  of  radar  site  above  equatorial  plane  is 
^  ““®CR  *  %  *  «SR> 


where 


*  (*‘81^  co89^,j^)(l-e)  tanQjjj^ 


Site  Centered  Certesian 


Elllpaold  centered  cartesian  inertial  coordinates  converted  to  ellipsoid* 
fixed  rectangular  with  origin  at  radar  (Fig.  10)  are  given  by  the  following: 


*  X  cos 


,'1'  .  y  .ini^ 


■  -X  sln'^  +  y  cos'i^ 


I  ««  z 

c 
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TBACKER  COOF(DIIiATES  ON  REET3iENCE  BLUFSOID 


FIGURE  10 


ffaaor  c;ciov«moas 


Site  Centered  Cagteftlaa  {ContM) 


vhejr« 


'i'  -  *  fi.  t 


**  *  ‘ 


w  ■ 


V  •  y* 

«  Wjj  « 


where  6  -  devletion  fron  local  ▼erticel.  end  eloB  •  A 

R  w  w 


u  -  u  coa*ia  ♦  ▼ 

«  n 


■  -U  +  V  CO«|j|j^ 


W  ■  V 


Pii<{e  33 


where  ■  an^le  throu^^  which  the  u  -  v  plaae  Is  rotated  to  Hx  radar 


aalQuth  xwference 


•  X  coe'i^  +  f  Mi3^  ♦ 

•  -i  eln^  ♦  coi^  -Q  •  x 


i  -  4 
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Bedar  Coxnreralons 


apherlcal  Radar  (Cont'd) 


A  ■  sin 


u  /Ui 

— TTTf 
+  v^)  J 


>),  0<A<2« 


E  -  cos*^  (|),  0<E<a 


R  «  g-  (uu  +  W  +  ww) 


:  vu  -  u'^ 

^  ”  ri — 57 

(u  +  V  ) 


_2,  2  2^ 

R  (u  +  V  ) 


Doppler  Radar 

Radar  cross -shaped  doppler  rates  cooputed  from  ellipsoid-fixed  rectangular 
coordinates  (Fig.  13) 


.  u 

P  ’  R 


<u 


(6,  -  S  )  ,  I-  (6^2  -  6/)  .  I  .  ,3) 


*  "3 

*i  '3^ 

■'^l 

^  *  1  'l'  -  1  'l^]  ' 
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Hatlar  Conversions 


Doppler  Radar  (Cont'd) 


S  (j  [(»!  •  *3)  *  J  -  63“^)  *  I  (»i5  -  Sj^)] 

-3  [-  »3  *  I '3^  4  »3l  -X  [-  *  I  »X^  *  I 'xfl 

S  [“  [<“x  -  '3’  ^  I  <''x‘’  -  ‘3")  *  I  (»x^  ■  »3^)] 


«  ■  ?  [“  [<*2  -  \>  *  I  <'a*  -  *  I  -  Bj,^)] 


“4 \ 4 B,3]  -»,[x*.,.a»^ 


2  +  1  ft  3 
2  ®2 


*  I  [('2  -  '4)  *  I  (Bj®  -  *  I  (Bj3  -  B^3)j 


Page  37 


Radar  Conversioiui 


Soppier  Radar  (Coat'd) 


ff  r*  [">2  -  »4>  *  I  o 


2  ft  Z\  ^ 

2 


where 


[l  ^  *  I  ^  I  \3]  -v^  (l  *  4  »2'  ^  I  »2^ ! 


(uu^  +  +  WWj)  -  I  +  M^) 


Redar  L-shaped  doppler  poeltlon  and  rates  computed  from  ellipsoid -fixed 
coordinates  (Fig.  13)  are  given  by  the  following.  Canstaats  are  delineated  on 
Page  ^2. 

9  9  9  0 

R”  «  u”  +  v“  +  w" 


R  +  -  uK 


1  -’''s] 


2 

D  j,  w  -  uir 


P  *  -  “*9  - ’*L0  - '«u] 


*16  * 


+  2wK, 


T»  •  H 
"©  ^  "1 


Radftr  Coaversions 


• '  r«(.v 


Doppler  Radar  (Cout'd) 


^^^7*  2vKi8  ^  2wKi^ 


u{u  -  Kg^^)  ^^  y(v  -  Kg^)  -f  v(v  -  Kgg) 


^  vKi3  -  ^  ^o 


Look  Ancles 


•  ,  ^17  ^^16  ~  ^  ^o 


Missile  look  sngles  from  radar  site  are 

+  (y  -  y. )  + 


BlAl  s  cos 


■if** 

T"^ — i - z  - 

L|(x  -  +  (y  -  y_)  +  (z  -  z  ) 

i  r  *  r  j 


RIA2  1=  tan 


-1 


(x  -  Xj.)C^  4  {}•  -  y^)  *  (z  -  Zj.) 

I*  ■  V’^x^'I'y- V>ty* 
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SITE  -  FIXED  RECTANGULAR  TO  SPHERICAL  CONVERSION 


FIGURE  12 
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DOPPIPR  RATE  STATION  SETUP 


FIGURE  13 


liaUMiiW  1-1  ■nirr— . .  . 


4l 


looli:  tofdes  (Cont'd) 
where 

^r  '  ^er  +  ^r) 


^r  “  V  ^  ffj,) 


z 

r 


z 

er 


Pitch  and  yww  look  an£les 
If  RIAl  =  90°  or  270° 


lAY  =  90°  with  opposite  sign  of  RIA2 
If  |rLA2  I  <  90°J  up  -•=  -90° 

If  |rU2  [  >  90°;  UP  «  90° 

If  RIAl  /  90°  or  270° 


tem  KUl 


sipRUl 

cosRUl 


UP  =  tan'^ 
UY  =  tan’^ 


iAJ 


[-( C08RIA2)  ( teujRUl) 
[-(elnRU2)  ( taaRIAl) 


UP  =  UP 


UY  -  UY 


If  RIAl  >  90° 

If  RU2  -90° 

UP  =  j  Sff  j 
UY  =  UY  +  180° 


i 


t. 


Lea  (Cont'd) 


If  RIA2  >  -90° 

If  RU2  ^  0° 

LAP  -  -  |i]ff  I 
lAY  «  Hy  +  180° 

If  RLA2  >  0° 

If  RIA2  ^  90° 

LAP  =  -  |Op| 
LAY  »  lAY  -  180° 
If  RIA2  >  90° 

LAP  - \m\ 
lAI  =  lAY  -  100° 


-4:/: 

Page  42 


I 


FIGURE  14 


Inltlallxatlon 


launch  Aalinuth  Estimation  (Cont'd) 

-  220.37 
«  20.621 
Kg  =  0. 

KSTO  =  11137.6 
KSTl  -  199.Wto6 
KST2  *  -.8J+8166 


K2ST0  =  lUl.675 
K2ST1  =  .0062625 
K2ST2  »  0. 

Launch  Attit\ide 


^  To  achieve  launch  vertical  the  b  , b  system  is  oriented  such  that 
£  points  along  th^x  direption,  7}  along  the  z  direction,  and  £  in  the  negative 


y  direction.  The  b  3  system  is  rotated  through  longitude  of  launch,  astronomic 
latitude  of  launch,  aiid  the  si^lement  of  launch  azimuth.  (Fig.  15) 


£  =  £  cos(-^  -  Q_t)  +  £  sin(-i^-  -  Q  t) 

Ij  C  ii  6 

n 

£  =  -  £  8ln(-^^  -  Q  t)  +  £  C08(-^-  -  Q  t) 

Li  C  C 


Inlti«Us«tloa 


XaSufiCu  Aw'trjL^Uu® 


tb«n 


i  m  i  co»(-«^)  -7}  •an(-o^) 

7}  m  b  »in(-0^)  *71  cso»(-Oj^) 


4’.r 

finaUy  ^ 

c  r 

7(^1  71  C08  (AZj^  -  Iflo)  -  C  •ln(AZj^  -  IBO) 
sin(itt;j^  -  100)  +  C  co»(AZ^  -  iflO) 

Launch  coordinate  rotational  quantity,  Hj^  la  eatlMated^  If  not  spaclfled^  aa 

Poaltlon  and  Velocity  Cociponenta  at  Launch 
Sea  level  ra&lua  of  launch  (Flg<  l6)  la 

1/2 


1  +  (1-e)^  ] 

1  +  (1-e)®  tan^«^  J 
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InitlallBAtion 


Poaition  and  Velocity  CoiaponentB  at  Launch  (Cont'd) 


Counponentfi  of  r  and  V  are 
Q  O 


*0  ■  '”•<^1,  *  [’■31^  “>"*CL  ’  <"1.  *  '"“dl] 

*>  [  'SI^  '•'■'•ct,  ♦  '•‘l  *  %l’  '“V] 

•“’cL  *  (4,  *  «8L>  '““dI, 


X  ®  -y 
o  ''o  e 


y_  =  X  Q 
o  o  e 


i  «  n 
~o 


The  siaiaLation  may  he  started  at  an  arbitrary  point  in  space  vith 


£  .  and  T} 


specified. 


If  the  position  and  velocity  components  are  given  in  a  reidar  coordinate 
system^  the  following  conversions  are  made; 


':Ei'  ' 


Pa*.**-  5-> 


•’t.^  <^y> 


Radar  to  CartesiBU  Coordinate  Conversion 


Ellipsoid  -  Spherical  Coordinates  and  Doppler  Rates  to  Fixed  Rectangular 


u  =  R  sinE  slnA 


V  =  R  sinE  cosA 


V  =  R  cosE 


r  2  1 

»„  -  ^  -  uK^  .  .Kg  -  vKjJ 

r  2  1 

=  R  +  Kq  -  uK^  -  vK^  -  wK^J 

Sg  =  ^R  +  K^2  "  ■  ''%lj 


vK,^  -  wK, 


^16  2(uK^^  +  v%4  ■«•  vKj^^) 


«o"«l 


Kgo  +  2(uKj^„  +  vKj^g  ->■  wK^q) 


K  +  no 

V/  \  ■  \4 


R  =  R 
o 


Rq  P  “  P  +  p(Rq  -  p) 


R^  4  =  R^  4  +  q(R^  -  4) 


D  .  KggCu-Kgi^)  KggCv-Kg^)  +  K2i(v-K26) 


IN 


P«ce 


wm. 


■iP" 
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Cuordinate  Conversion  (Coat'd) 


[ki5(v.K25)-K,^(v^6>]  "  ^  -T^  ' 


*  *  L  [''l3^""'^6^"'^5^“"'^4^]  *  T5“[*19^'*"*24^"''i7^''"*26^]  D  ’  *22 


*  '  ¥  [*i4(“-*21.)-*13<’'-*25’]  ♦  ¥[*17<''-*25’-*i8<“-*24>]  *  ■  *21 


Dealer  radar  coefficients  (see  Fig.  I3). 


2(Vo) 


2(Vo) 


(u^)^  +  (Vq)^ 


K5  -  2{u3^) 

Kg  =  2(vj^) 

-  2(w^) 


(Uj^)^  +  (v^)^  +  (Wj^)® 
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Inltiftiiaation 


Radar  to  Cartes ian  Coordinate  ConvorBlon  (Cont'd) 


=  2(^2) 

Kii  = 

''12  •  (“2)"  -  (’'2)' 


(V2)"  "  ^2' 


'‘13  =  “1  -  “o 

■'ll.  =  ''1  -  ''o 

''15  ■  '*1  •  ”0 

V  H  2  H  2 

^16  ~ 


•^17  =  % 

^18  =  ^2  - 
"'ig  '  ''2  - 

»r  2  .  2 

*^20  “  ■  ^2 


^21  =  ^13’'i8  -  Vlh 
^22  "  ^15*^17  ■  ^19^13 
^23  "  “  ^18^15 


=  “o 
^25  == 


^26  "  ''o 


Inltiall  utioD 


Radar  to  Cartesian  Coordljmte  Conversion  (Cont'd) 


0  cobO, 


C08(Q^ 


r  ■  j6in(Q^  ^ 


-BlnCQ^  t  + 

cob(Q^  t  +  6^j^) 


-I  L. 


“  T  p  •!  p 

co8(Q^  t  +  -ain(W^  t  +  0  S  \-^ 

ir  -  aln(n^  t  +  coo(Qj  t  +  B^jj)  0  x 


L°J 


Attitude  from  Angles  of  Attack 

The  following  conversions  are  made  to  Initialize  the  attitude  of  the  missile 
when  OL  and  a  are  Imovn: 


Inltl&llz»tloD 


Attitude  fi'om  Angles  of  Attack  (Cont’d) 


CMP " f *  •  f  '  i'  * '  ■  £  1 

I  cl  I  % 


8inP  =  1  -  cos  3 


^  =  ~  -  cosp  •  £ 

o  SqJ  slop 


>'  >'  ' 

i 

o  'o 


C  "  [  ®°®“p  ■  £  ■  ^  1  ’  £ 

^  o  O'*  o 

■  £^  "■  \ 
i-X-  \ 


Pocc  5^' 


Initifilitation 

Tie  Do\m 

Derivatives  of  equotions  of  motion  to  maintain  vertical  attitude  during 
launch  hold  down  are  coinputed  from 

•  n, 

“t  •  -  ( ^i)  • 

“V  ■  ■ 
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APPENDIX  A 
Symbol  Definitions 

Semi -major  axis  of  reference  ellipBoid. 

Azimuth  of  launch,  measured  clockwise  from  north. 

Wind  azimuth,  measured  clockwise  fi-om  north. 

Semi -minor  axis  of  reference  ellipsoid. 

Speed  of  Eoiuid. 

'j'abular  dree  coefficient. 

A  ^ 

Ellipticity  of  reference  ellipsoid  =  . 

A 

C  — 

Eccentricity  =  — j  where  C  is  the  distance  from  the  center 

A 

of  the  ellipsoid  to  a  focus  along  the  major  axis. 

Gravitational  Acceleration. 

Longitude. 

Astronomic  latitude,  angle  between  the  equitorial  plane  and 
the  local  direction  of  gravity. 

Geocentric  latitude,  angle  between  the  equatorial  plane  and  . 

Geodetic  latitude,  angle  between  the  eq\iatorial  plane  and  the 
normal  to  the  reference  ellipsoid. 

Multiplier  for  tabular  thrust  F^. 

VJlnd  perturbation  multiplier. 

Multiplier  for  t^ular  maae  flow  rate  M^. 

Multiplier  for  tabular  pitch  rate. 

Multiplier  for  tabular  roll  rate. 

Multiplier  for  tabular  yaw  rate. 
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>iPPENDIX  A  (Cont'cl) 

Macs  01'  miscile. 

Refei'cncc  Initial  mass  for  ablation  corqputati on . 

Mach  number. 

Atmospheric  pressure. 

Dynamic  pressure. 

Atmospheric  density. 

Missile  position  vector  whose  components  are  x,y,z. 
Average  radius  of  geoid  (spherical  approximation). 

Missile  cross-sectional  area. 

Atmospheric  temperature. 

Longitude . 

Tabular  wind  velocity. 

Components  of  missile's  velocity  vectoi'  V. 

Reference  ellipsoid's  rotation  rate. 

Angular  rotational  rate  about  respectively  pitch,  yav  and 
roll  axes. 

Respectively  pitch,  yaw  end  roll  rates  from  the  guidance 
program . 

Respectively  pitch,  yaw  and  roll  drifts. 
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Two  general  schemes  of  iteration  are  used.  The  first,  called  (for  lack  of  a 

better  name)  "paired  function",  requires  that  each  parameter  to  be  determined 
be  paired  with  a  function  to  be  constrained.  The  second  form  of  iteration, 
termed  "matrix",  does  not  require  pairing  of  the  functions  and  parameters  but 
docs  not  allow  minimizing  or  maximizing  of  functions. 

In  either  type  of  iteration,  control  of  the  portions  of  a  trajectory  over  vdiich 
iteration  is  to  take  place  is  provided.  The  control  system  allows  "nesting"  of 
"paired  function"  sets  within  other  "paired  function"  or  "matrix"  sets.  This 
capability  allows  combination  of  the  strong  points  of  the  two  general  schemes 
according  to  the  needs  of  the  individual  problem. 

Paired  Functions 

After  an  initial  t-^“;*«ctory  is  integrated,  a  perturbation  is  added  to  the  para¬ 
meter  in  questioi.  end  t'other  trajectory  is  integrated.  Inverse  La  Granglan 
interpolation  or  extr  is  then  used  to  estimate  a  new  valye  for  the 

parameter.  (The  desii-'d  value  of  the  function  is  used  as  argument.)  The  traj¬ 
ectory  is  again  integrated  and  the  results  enable  the  use  of  a  higher  order 
polynomial  in  tlie  Interpolation  or  extrapolation  procedure.  The  process  is  used 
repeatedly  until  convergence.  The  polynomial  order  may  be  limited  by  input,  or 
in  any  case  at  order  5. 

If  minimization  or  maximization  is  called  for,  a  second  perturbation  is  added  to 
the  parameter  and  a  third  trajectory  is  produced.  La  Grangian  interpolation  is 
then  used  in  a  "step -and -bisect"  iterative  procedure  \diich  determines  minimum  or 
maximum  for  the  polynomial  fitting  the  available  points.  If  no  minimum  or  maxi¬ 
mum  is  available,  another  perturbation  is  added  to  the  parameter.  In  either  case 
another  trajectory  is  produced  and  the  new  point  is  used  to  enable  a  higher  order 
polynomial  to  be  used  in  the  "step-end -bisect"  solution.  Again  the  La  Grangian 
polynomial  order  may  be  limited  by  input,  or  at  order  5- 
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I  Paired  Functions  (cont'd) 

The  'btep>and'-hi8ect"  process  is  said  to  converge  vhen  the  improvement  in  the 
function  is  less  than  one  half  the  (input)  accuracy  requirement.  The  function 
1  maximization  (minimization)  is  complete  vhen  the  difference  "between  output  of 
I  the  "step-and-hisect"  procedure  and  the  previous  value  of  the  function  is  less 
then  the  (input)  accuracy  requirement. 


MATRIX  Iteration 

The  Newton -Raphson  method  for  iterative  linear  differential  correction  may  be 
stated  in  its  matrix  form: 


(1) 


Where  Ay  !  consists  of  the  differences  fr<xn  desired  conditions,  andjAxJ 
is  the  GUGgus-ucd  chcnccc  to  the  pcremeters  in  question. 

The  icvtisl  uorivotive  elements  ^1  may  be  replaced  by  their  finite  difference 
cqulvcJ.ents ,  A  1,J  ,  In  the  N-CiAG^  pri  'rtni,  in  order  to  enable  the  use  of  the 


Genercllzed  Secant  Method,  the  resultant  equation  is  re-written  as  two  equations  in 
the  <-i'ui>folcnt  form: 
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AX,  0.  .  -  .  . 

(3) 

Ax 

m 

± 

? 

0. 

• 

a 

*  A  X 

•  u.* 

! 

—  -j 

Where  A  ^ j  indicates  the  difference  in  produced  by  a  change  A  in  X j 

Hie  elements  AX^  a  we  obtained  by  repeatedly  Integrating  the  trajectory  with  c 
each  of  the  parameters  in  question  perturbed  in  its  turn.  The  elements  of 
the  diagonal  matrix  (3),  are  the  perturbation  values.  The  elements  ofjAT^end 
|Ay^  ^  are  divided  by  normalizing  elements  (which  are  input  quantitie^  in 
order  to  scale,  or  establish  dimensional  homogeneity. 

The  first  equation  is  solved  for  J  in  a  "least  squares"*  sense,  which  ^lovs 
the  general  case,  m/n.  The  multiplication  implied  in  the  second  equation  is  then 
perfc^d, yielding  the  eet  [a0  ,  the  elet.ents  of  .hlch  ere  added  to  the  appro- 
priate  parameters. 


After  a  solution  of  the  above  sort  is  accomplished,  and  the  trajectory  representing 
the  suggested  parameters  integrated  and  found  to  be  still  unconverged,  the  Generalized 
Secant  Method  may  be  utilized  to  improve  the  convergence  without  integrating  the  num¬ 
ber  of  trajectories  (n),  necessary  to  recompute  the  matrices  for  equations  (2)  and 
(3).  In  this  method  a  new  column,  consisting  of  the  negative  of  the  normalized 
(scaled)  deviations  from  conditions  at  the  previous  nominal  trajectory  integration, 
is  appended  to  each  m£'.:rix  as  the  first  column,  and  is  added  to  each  of  the  remaining 
columns,  producing  matrices  with  one  more  column  than  the  previous  matrices. 


-Ay 

■A'^2,0 

;  M 

1 

.  .  •  A  Yjn^P  'AYui^q 

1 

» 

L.  _ 

-  i-wiijI+F,  Ger.crul  Linetr  Equation  Solver  by  D.  D.  Morrison,  Sept.  2h,  I959  (CDRC 

■  rr’.pit.rr  r>ub:' ,  u t Ine )  . 
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(5) 


Ax 


« 


•  • 


-AX 


n,o 


The  procesn  ma>  be  used  repeatedly,  each  time  defining  the  previous  trajectory  to 
lie  Lhcj  vrevlous  nominal  trajectory. 

The  last  column  may  be  deleted  after  any  number  of  additional  colvmns  have  been 
accumulati'id. 


A  satisfactory  rationale  for  deletion  or  non-deletion  of  the  last  column  is  not 
presently  available,  but  experience  indicates  that  the  last  column  should  be  kept 
if  the  number  of  functions  is  approximately  equal  to  the  number  of  parameters  (msn), 
and  deleted  vhen  the  number  of  functions  is  markedly  greater  than  the  number  of 
pcrcineters  (n>n).  Tests  vith  m  and  n  initially  equal  to  3  have  shovn  a  20^6  Increase 
in  elflciency  ^en  the  matrices  are  allowed  to  expand  indefinitely. 


In  particular,  when  Inaccuracies  due  to  system  non-linearities  (roiand-off  error 
accumulations,  integration  truncation  error,  staging  condition  annomalies,  etc.) 
oi^e  on  ohc  order  of  the  iteration  convergence  accuracies  required,  it  has  been  found 
necessrry  to  l Ilow  the  matrices  to  expand  in  order  to  gain  convergence. 


An  addii^ionol  advantage  of  the  Generalized  Secant  Method  is  that  perturbation  delta 
magnitudes  a  less  dominant  role  in  convergence  efficiency.  In  many  cases 

converger.eo  s  worsened  by  a  Newton -Raphson  solution  because  of  erroneous  partial 
derivative  Information  due  to  improper  perturbation  sizes.  In  such  cases  it  is  not 
uncommon  to  find  that  the  Generalized  Secant  Method  can  proceed  to  convergence. 


The  Generalized  Secant  Method  is  used  repeatedly,  until  the  gain  in  one  step  is  less 
ihun  Input  quantity  SCGN.  (The  gain  is  measured,  as  the  ratio  of  the  summations 


squares  of  the  elements  cf  the  normalized  error  vecto: 


{Af]  , 


for  two  successive 


.'aiu 


atlon  passes.)  Thus,  if  any  solution  pass  (whether  the  solution  comes 


Best  Available  Copy 
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from  Nevton^apbson  or  Oeneralised  Secant  Kothod)  la  not  materially  better  than  the 
previous  solution  or  ncnlnal  passj  the  progrsai  reecmputes  the  matrices  by  perturba¬ 
tion  methods,  ffote.  that  If  SOfflf  is  a  very  large  maibw,  the  Oeneralised  Secant 
Method  Is  not  attempted. 

Eigenvalues  associated  with  the  solution  of  eqmation  (2)  or  (U)  are  printed  foUovlng 
each  evaluotlon.  These  Eigenvalues  are  of  use  only  if  it  Is  understood  that  they-  are 
for  the  matrix  transpose  multiplied  by  the  matrix^  and  it  Should  be  remembered  that  the 
elements  are  scaled.  The  ratio  of  the  largest  Eigenvalue  to  the  smallest  (knovn  as  the 
condition  number)  is  sometimes  a  clue  to  the  degree  of  linear  dependence  in  the  matrix. 

A  possibly  confusing  issue  is  that  the  overall  perturbation  in  the  column  associated 
"ith  a  particular  Eigenvalue  also  affects  the  condition  number.  Zero  Eigenvalues 
should  not  be  found  for  a  Nevton-Raphson  solution,  but  are  often  found  for  a  Generalized 
Secant  Method  solution.  In  fact,  usually  only  (n)  Eigenvalues  are  large  enough  to  have 
any  effect,  even  thoucdi  the  matrix  has  been  alloved  to  expand  considerably  larger  than 
the  original  (n)  columns. 

In  any  cose,  the  Eigenvalues  are  associated  vith  input  parameters  in  reverse  order  to 
their  printing.  That  is,  the  last  Eigenvalue  is  associated  with  the  first  iteration 
variable. 

^  matrix  iteration  is  determined  to  have  converged  on  either  of  tvo  criteria. 

1.  The  error  in  each  of  the  functions  is  less  than  the  corresponding  accuracy 
requirement . 

2.  The  change  In  the  error  (for  each  of  the  functions  ^Ich  fails  the  first 
criterion)  is  less  than  one  half  the  corresponding  accuracy  retjuirement. 

The  change  is  measured  from  the  previous  solution  or  the  nominal  trajectory. 
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